##### Step 1: load data and packages #####



library(dplyr)
library(magrittr)
library(readr)
library(stringr)
library(ggplot2)
library(fixest)
library(stargazer)
library(doMC)
library(haven)
library(data.table)

options(scipen = 5)

base_dir <- "..." # Manually change home folder 
data_dir <- paste0(base_dir, "data/")

fixest_form <- function(
  y,
  x, 
  fe=NULL,
  x_re = FALSE,
  df = NULL, 
  endog = NULL,
  insts = NULL
) {
  
  if (x_re & !is.null(df)) {
    x <- x %>% map(~ grep(., colnames(df), value=T)) %>% reduce(c)
  }
  
  if (!is.null(endog)) {
    ivleft <- paste(endog, collapse = " + ")
    ivright <- paste(insts, collapse = " + ")
    ivpart <- paste(ivleft, ivright, sep = " ~ ")
  }
  
  xpart <- paste(x, collapse = " + ")
  if (!is.null(fe)) {
    fepart <- paste(fe, collapse = " + ")
  } 
  
  if(is.null(fe) & is.null(endog)) {
    y %>% 
      paste(xpart, sep=" ~ ") %>%
      as.formula
  } else if (is.null(fe)) {
    y %>% 
      paste(xpart, sep=" ~ ") %>%
      paste(ivpart, sep = " | ") %>% 
      as.formula
  } else if (is.null(endog)) {
    y %>% 
      paste(paste(xpart, fepart, sep="|"), sep=" ~ ") %>%
      as.formula
  } else {
    y %>% 
      paste(paste(xpart, fepart, sep="|"), sep=" ~ ") %>% 
      paste(ivpart, sep = " | ") %>% 
      as.formula
  }
  
  
}

tablestyle  = style.tex(main="aer",
                        fixef.suffix = " FEs",
                        fixef.where ="var",
                        fixef.title = "",
                        stats.title = "\\midrule",
                        tablefoot=F,
                        yesNo="\\checkmark")

setwd(paste0(data_dir, "/gfk"))

files <- list.files()
files <- files[str_sub(files, 11, 11)=="0"]
print(files)
csv.files <- files[str_sub(files, -3, -1)=="csv"]

for(file in csv.files) {
  assign(str_sub(file, 1, -5), read.csv(file))
}

for(file in str_sub(files, 1, -5)) {
  assign(paste0("dict_", file), names(get(file)))
}

# Interview dates
interview_dates <- read.csv("mri_respid_zip_date.csv")

# Channel positions
chpos <- readRDS(paste0(data_dir, "cable_zip.rds"))

# zip demog
load(paste0(data_dir, "zip_demographics_2000.RData"))





##### Step 2: combine data within survey wave #####



mediamark_01 <- fnctp_gfk_01

mediamark_03 <- fnctp_gfk_03

merge.vars.05 <- intersect(dict_fnctp_gfk_05, dict_fnctp_gfk_05_cable)
mediamark_05 <- left_join(fnctp_gfk_05, fnctp_gfk_05_cable, by=merge.vars.05)
stopifnot(nrow(mediamark_05)==nrow(fnctp_gfk_05))

merge.vars.07 <- intersect(dict_fnctp_gfk_07, dict_fnctp_gfk_07_politics)
mediamark_07 <- left_join(fnctp_gfk_07, fnctp_gfk_07_politics, by=merge.vars.07)
stopifnot(nrow(mediamark_07)==nrow(fnctp_gfk_07))

merge.vars.07.cable <- intersect(names(mediamark_07), dict_fnctp_gfk_07_cable)
mediamark_07 <- 
  left_join(mediamark_07, fnctp_gfk_07_cable, by=merge.vars.07.cable)
stopifnot(nrow(mediamark_07)==nrow(fnctp_gfk_07))

merge.vars.09 <- intersect(dict_fnctp_gfk_09, dict_fnctp_gfk_09_cable)
mediamark_09 <- left_join(fnctp_gfk_09, fnctp_gfk_09_cable, by=merge.vars.09)
stopifnot(nrow(mediamark_09)==nrow(fnctp_gfk_09))

years <- ls(pattern="mediamark_")
years <- str_sub(years,-2,-1)





##### Step 4: merge in zip codes and survey years by RespID #####



mediamark_01 %<>% 
  left_join(interview_dates[, c("RespID", "zipcode", "int_year")])
mediamark_03 %<>% 
  left_join(interview_dates[, c("RespID", "zipcode", "int_year")])
mediamark_05 %<>% 
  left_join(interview_dates[, c("RespID", "zipcode", "int_year")])
mediamark_07 %<>% 
  left_join(interview_dates[, c("RespID", "zipcode", "int_year")])
mediamark_09 %<>% 
  left_join(interview_dates[, c("RespID", "zipcode", "int_year")])





##### Step 5: merge in channel position by zip code #####



chpos %<>% mutate(zipcode = as.integer(zipcode))
mediamark_01 %<>% left_join(chpos, by=c("zipcode", "int_year"="position_year"))
mediamark_03 %<>% left_join(chpos, by=c("zipcode", "int_year"="position_year"))
mediamark_05 %<>% left_join(chpos, by=c("zipcode", "int_year"="position_year"))
mediamark_07 %<>% left_join(chpos, by=c("zipcode", "int_year"="position_year"))
mediamark_09 %<>% left_join(chpos, by=c("zipcode", "int_year"="position_year"))





##### Step 5: filter to cable-only viewers #####



for(year in years) {
  
  dat <- get(paste0("mediamark_", year))
  
  cable_pos = "X.13015"
  satellite_pos = "X.13029"
  
  cable.var <- names(dat)[which(str_detect(names(dat), cable_pos))]
  print(cable.var)
  
  satellite.var <- names(dat)[which(str_detect(names(dat), satellite_pos))]
  print(satellite.var)

  rows.to.keep <- which(dat[[cable.var]]==1 & dat[[satellite.var]]==0)
  dat <- dat[rows.to.keep, ]
  
  assign(paste0("mediamark_", year), dat)
  
}





##### Step 6: create political activity measures #####



for(year in years) {
  
  dat <- get(paste0("mediamark_", year))
  
  if(year=="01") {
    activity_start = "X.100511"
    activity_end = "X.100523"
    voting = "X.100511"
    fundraising = "X.100510"
    rally = NA_character_
  }
  if(year=="03") {
    activity_start = "X.126341"
    activity_end = "X.126353"
    voting = "X.126341"
    fundraising = "X.126340"
    rally = NA_character_
  }
  if(year=="05") {
    activity_start = "X.09441"
    activity_end = "X.09455"
    voting = "X.09441"
    fundraising = "X.09454"
    rally = NA_character_
  }
  if(year=="07") {
    activity_start = "X.09441"
    activity_end = "X.09455"
    voting = "X.09441"
    fundraising = "X.09454"
    rally = "X.09446"
  }
  if(year=="09") {
    activity_start = "X.09441"
    activity_end = "X.09455"
    voting = "X.09441"
    fundraising = "X.09454"
    rally = "X.09446"
  }
    
    activity.vars <- 
      names(dat)[
        which(str_detect(names(dat), activity_start)):
          which(str_detect(names(dat), activity_end))]
    print(activity.vars)
    
    activity_aggregate <- paste0(
        "dat %<>% rowwise() %>% mutate(
    any_political_activity = max(",
        paste0(activity.vars, collapse=", "),
        "  ))")
    cat(activity_aggregate)
    eval(parse(text=activity_aggregate))
    
    dat$voting <- dat[[which(str_detect(names(dat), voting))]]
    dat$fundraising <- dat[[which(str_detect(names(dat), fundraising))]]
    
    if(!is.na(rally)) {
      dat$rally <- dat[[which(str_detect(names(dat), rally))]]
    }
    if(is.na(rally)) {
      dat$rally <- NA_integer_
    }
    
  assign(paste0("mediamark_", year), dat)
  
}





##### Step 7: measure ideology #####



for(year in years) {
  
  dat <- get(paste0("mediamark_", year))
  
  if(year=="01") {
    ideology_start = "X.99551"
    ideology_end = "X.99555"
  }
  if(year=="03") {
    ideology_start = "X.125381"
    ideology_end = "X.125385"
  }
  if(year=="05") {
    ideology_start = "X.139131"
    ideology_end = "X.139135"
  }
  if(year=="07") {
    ideology_start = "X.149651"
    ideology_end = "X.149655"
  }
  if(year=="09") {
    ideology_start = "X.157681"
    ideology_end = "X.157685"
  }
  
  ideology.vars <- 
    names(dat)[
      which(str_detect(names(dat), ideology_start)):
        which(str_detect(names(dat), ideology_end))]
  print(ideology.vars)
  
  ideology.dat <-
    apply(dat[,ideology.vars], 2, sum)
  ideology.dat %<>% as.data.frame() 
  names(ideology.dat) <- "count"
  
  ideology.dat$rownum <- 1:nrow(ideology.dat)
  ideology.dat %<>% mutate(
    ideology_quintile = -(rownum-mean(ideology.dat$rownum))
  )
  
  # 2022-06-07: change ideology to 3-tier -----
  ideology.dat %<>% mutate(
    ideology_quintile = case_when(
      ideology_quintile < 0 ~ -1,
      ideology_quintile == 0 ~ 0,
      ideology_quintile > 0 ~ 1,
      TRUE ~ NA_real_
    )
  )
  
  for(var in rownames(ideology.dat)) {
    ideology.dat[[var]] <- 
      as.integer(var==rownames(ideology.dat))
  }
  
  dat %<>% 
    left_join(
      ideology.dat[, c(rownames(ideology.dat), "ideology_quintile")], 
      by=rownames(ideology.dat))
  
  assign(paste0("mediamark_", year), dat)
  
}





##### Step 8: condense individual demographic variables #####




### HH income: take midpoint of income bracket ###

# HH income: X.08191 -- X.08205

for (year in years) {
  
  dat <- get(paste0("mediamark_", year))
  
  if(year=="01") {
    income_start = "X.08191"
    income_end = "X.08203"
  }
  if(year %in% c("03", "05", "07", "09")) {
    income_start = "X.08191"
    income_end = "X.08205"
  }
  
  income.vars <- 
    names(dat)[
      which(str_detect(names(dat), income_start)):
        which(str_detect(names(dat), income_end))]
  print(income.vars)
  
  income <- dat[,income.vars] %>% apply(2, sum)
  income %<>% as.data.frame() 
  names(income) <- "count"
  
  temp <- strsplit(rownames(income), "\\.")
  bin_start <- sapply(temp, function(x) {
    return(as.integer(x[length(x)-2]))
  })
  bin_end <- sapply(temp, function(x) {
    return(as.integer(x[length(x)-1]))
  })
  
  income %<>% bind_cols(bin_start) %>% bind_cols(bin_end)
  names(income) <- c("count", "bin_start", "bin_end")
  
  income %<>% mutate(
    bin_end = ifelse(is.na(bin_end), Inf, bin_end)
  )
  
  avglnparms <- c(3.615872,0.798343) # parameters for the average distribution (32% / 30% / 28% / 8% / 2%)
  
  lognorm.incvals <- function(par, boundaries) {
    mu <- par[1]
    sigma <- par[2]
    expvals <- apply(boundaries, 1, function(bounds){
      expval <- integrate(function(x) {x * dlnorm(x, mu, sigma)}, bounds[1], bounds[2])	
      expval$value / (plnorm(bounds[2], mu, sigma) - plnorm(bounds[1], mu, sigma))
    })
    expvals
  }
  
  income$hh_income <- 
    lognorm.incvals(
      avglnparms, 
      boundaries=as.matrix(income[, c("bin_start", "bin_end")]/1000))*1000
  
  for(var in rownames(income)) {
    income[[var]] <- as.integer(var==rownames(income))
  }
  
  dat %<>% left_join(income[, c(rownames(income), "hh_income")], by=rownames(income))
  
  assign(paste0("mediamark_", year), dat)

}



### Respondent age ###

# *Respondent* age: X.04311 -- X.04324 

for(year in years) {
  
  dat <- get(paste0("mediamark_", year))
  
    age_start = "X.04311"
    age_end = "X.04324"

    age.vars <- 
      names(dat)[
        which(str_detect(names(dat), age_start)):
          which(str_detect(names(dat), age_end))]
    print(age.vars)
    
    age.dat <-
      apply(dat[,age.vars], 2, sum)
    age.dat %<>% as.data.frame() 
    names(age.dat) <- "count"
    
    temp <- strsplit(rownames(age.dat), "\\.")
    bin_start <- sapply(temp, function(x) {
      return(as.integer(x[length(x)-2]))
    })
    bin_end <- sapply(temp, function(x) {
      return(as.integer(x[length(x)-1]))
    })
    age.dat <- 
      bind_cols(
        age.dat,
        bin_start,
        bin_end
      )
    names(age.dat) <- 
      c("count", "bin_start", "bin_end")
    
    age.dat %<>% mutate(
      cum_freq = cumsum(count)/sum(count)
    )
    
    age.dat %<>% mutate(
      age_quintile = case_when(
        cum_freq <= 0.5 ~ 0,
        cum_freq > 0.5 ~ 1,
        TRUE ~ NA_real_
      )
    )
    stopifnot(sum(is.na(age.dat$cum_freq))==0)
    
    for(var in rownames(age.dat)) {
      age.dat[[var]] <- 
        as.integer(var==rownames(age.dat))
    }
    
    dat %<>% 
      left_join(
        age.dat[, c(rownames(age.dat), "age_quintile")], 
        by=rownames(age.dat))
  
    dat %<>% mutate(
      age_10s = 
        (`X.04311...Respondent...Age..18.Years`==1 | 
           `X.04312...19.Years`==1),
      age_20s =
        (`X.04313...20.Years`==1 |
           `X.04314...21.Years`==1 |
           `X.04315...22.24.Years`==1 |
           `X.04316...25.29.Years`==1),
      age_30s =
        (`X.04317...30.34.Years`== 1 |
           `X.04318...35.39.Years`==1),
      age_40s = 
        (`X.04319...40.44.Years` == 1 |
           `X.04310...45.49.Years`==1),
      age_50s =
        (`X.0431x...50.54.Years` == 1 |
           `X.0431y...55.59.Years`==1),
      age_60s =
        (`X.04321...60.64.Years` == 1 |
           `X.04322...65.69.Years`==1),
      age_70s_80s =
        (`X.04323...70.74.Years` == 1 |
           `X.04324...75..Years`==1)
    )
    
  assign(paste0("mediamark_", year), dat)

}




### Black/White/Hispanic ###

for(year in years) {
  
  dat <- get(paste0("mediamark_", year))

  if(year %in% c("01", "03")) {
    white_pos = "X.08361"
    black_pos = "X.08362"
    hispanic_pos = "X.08401"
  }
  if(year %in% c("05", "07", "09")) {
    white_pos = "X.07071"
    black_pos = "X.07072"
    hispanic_pos = "X.08401"
  }
  
  white.var <- names(dat)[str_detect(names(dat), white_pos)]
  black.var <- names(dat)[str_detect(names(dat), black_pos)]
  hispanic.var <- names(dat)[str_detect(names(dat), hispanic_pos)]
  print(white.var)
  print(black.var)
  print(hispanic.var)
  
  dat$white <- dat[[white.var]]
  dat$black <- dat[[black.var]]
  dat$hispanic <- dat[[hispanic.var]]
  
  assign(paste0("mediamark_", year), dat)
  
}



### Respondent education ###

# *Respondent* education: X.05751 -- X.05757

for(year in years) {
  
  dat <- get(paste0("mediamark_", year))
  
  grad_pos = "X.05755"
  postgrad_pos = "X.05756"
  
  grad.var <- names(dat)[str_detect(names(dat), grad_pos)]
  postgrad.var <- names(dat)[str_detect(names(dat), postgrad_pos)]
  print(grad.var)
  print(postgrad.var)
  
  dat$college_degree <- as.integer(dat[[grad.var]]==1 | dat[[postgrad.var]]==1)
  
  assign(paste0("mediamark_", year), dat)
  
}



### HOH gender ###

for(year in years) {
  
  dat <- get(paste0("mediamark_", year))
  
  male_pos = "X.04481"
  
  male.var <- names(dat)[str_detect(names(dat), male_pos)]
  print(male.var)
  
  dat$man <- dat[[male.var]]
  
  assign(paste0("mediamark_", year), dat)
  
}



### Viewership ###

for(year in years) {
  
  dat <- get(paste0("mediamark_", year))

  for(channel in c("cnn", "fnc", "msnbc")) {
  
    if(channel=="cnn") {
      hour_start = "X..13341"
      hour_end = "X..13340"
      any = "X..1334x"
    }
    if(channel=="fnc") {
      hour_start = "X..13561"
      hour_end = "X..13560"
      any = "X..1356x"
    }
    if(channel=="msnbc") {
      hour_start = "X..13781"
      hour_end = "X..13780"
      any = "X..1378x"
    }
    
    channel.vars <- 
      names(dat)[
             which(str_detect(names(dat), hour_start)):
               which(str_detect(names(dat), hour_end))]
    print(channel.vars)
    any.var <- 
      names(dat)[which(str_detect(names(dat), any))]
    print(any.var)
    
    channel.dat <-
           apply(dat[,channel.vars], 2, sum)
    channel.dat %<>% as.data.frame() 
    names(channel.dat) <- "count"
    
    temp <- strsplit(rownames(channel.dat), "\\.")
    bin_start <- sapply(temp, function(x) {
      return(as.integer(x[length(x)-3]))
    })
    bin_end <- sapply(temp, function(x) {
      return(as.integer(x[length(x)-2]))
    })
    channel.dat <- 
      bind_cols(
        channel.dat,
             bin_start,
             bin_end
           )
    names(channel.dat) <- 
      c("count", "bin_start", "bin_end")
    
    command <- paste0(
      "channel.dat %<>% mutate(",
      paste0(channel,"_hours"), " = case_when(
          is.na(bin_start) ~ as.numeric(bin_end),
          is.na(bin_end) & bin_start < 25 ~ 25,
          !is.na(bin_start) & !is.na(bin_end) ~ (bin_start+bin_end)/2,
          TRUE ~ NA_real_
        )
      )"
    )
    eval(parse(text=command))
    
    for(var in rownames(channel.dat)) {
      channel.dat[[var]] <- 
        as.integer(var==rownames(channel.dat))
    }
    
    dat %<>% 
      left_join(
        channel.dat[, c(rownames(channel.dat), paste0(channel,"_hours"))], 
        by=rownames(channel.dat))
    stopifnot(!(1 %in% dat[is.na(dat[[paste0(channel,"_hours")]]), rownames(channel.dat)]))
    dat[[paste0(channel,"_hours")]] <- 
      ifelse(is.na(dat[[paste0(channel,"_hours")]]), 0, dat[[paste0(channel,"_hours")]])
    
    dat[, paste0(channel, "_any")] <- dat[[any.var]]
  
  }
  
  assign(paste0("mediamark_", year), dat)
  
}





##### Step 9: Combine mediamark data across survey waves #####



vars.to.keep <- 
  c("RespID", "zipcode", "int_year", "has_MSNBC", "has_Fox.News.Channel", 
    "chan_config", "pos_lin_fnc", "pos_lin_msnbc", "num_channels", "num_bc", 
    "subs", "stid",
    "any_political_activity", "voting", "fundraising", "rally",
    "ideology_quintile",
    names(mediamark_01)[str_detect(names(mediamark_01), "age_")],
    "hh_income", "white", "black", "hispanic", "college_degree", "man", 
    "cnn_hours", "cnn_any", "fnc_hours", "fnc_any", "msnbc_hours", "msnbc_any")

mediamark <- bind_rows(
  mediamark_01[, vars.to.keep],
  mediamark_03[, vars.to.keep],
  mediamark_05[, vars.to.keep],
  mediamark_07[, vars.to.keep], 
  mediamark_09[, vars.to.keep])





##### Step 10: merge in zip code variables #####



zip_controls <- c(
  "pct_black", "pct_asian", "pct_other", "pct_hisp",
  "pct_male",
  "pct_age_10s", "pct_age_20s", "pct_age_30s", "pct_age_40s", "pct_age_50s", "pct_age_60s", "pct_age_70s", "pct_age_80s",
  "i(income_dec)",
  "pct_hsgrad", "pct_somecollege", "pct_bach", "pct_postgrad",
  "pct_suburban", "pct_urban",
  "totpop"
)

census <- census2000
head(census$zipcode)
census %<>% mutate(
  zipcode = as.integer(zipcode)
)
census %<>% group_by(zipcode) %>% slice(1) %>% ungroup()

# Construct income_dec_zip

income.vars <- 
  names(census)[
    which(names(census)=="inc_less_10k"):
      which(names(census)=="inc_200_plus")]

temp <- census %>% dplyr::select(zipcode, income.vars)
for(var in income.vars) {
  sum.vars <- income.vars[1:which(income.vars==var)]
  colsum <- paste0(
    "temp[['", paste0("cum_", var), "']] <- temp[['",
    paste0(sum.vars, collapse= "']] + temp[['"),
    "']]"
  )
  cat(colsum)
  eval(parse(text=colsum))
}

cum.vars <- names(temp)[str_sub(names(temp), 1, 4)=="cum_"]

temp$inc_med_cat <- NA_character_
for(i in 1:nrow(temp)) {
  if(is.na(temp$cum_inc_200_plus[i])) {
    next
  }
  if(!is.na(temp$cum_inc_200_plus[i])) {
    med.cat.pos <- which((temp[i, cum.vars]-0.5)^2==min((temp[i, cum.vars]-0.5)^2))
    med.cat = cum.vars[med.cat.pos]
    temp$inc_med_cat[i] <- med.cat
  }
}

temp$inc_med_cat = factor(
  temp$inc_med_cat,
  levels=c(
    "cum_inc_less_10k",
    "cum_inc_10_15",
    "cum_inc_15_20",
    "cum_inc_20_25",
    "cum_inc_25_30",
    "cum_inc_30_35",
    "cum_inc_35_40",
    "cum_inc_40_45",    
    "cum_inc_45_50",    
    "cum_inc_50_60",    
    "cum_inc_60_75",   
    "cum_inc_75_100",
    "cum_inc_100_125",
    "cum_inc_125_150",
    "cum_inc_150_200"
  )
  )

temptable <- table(temp$inc_med_cat) %>% prop.table()
temptable %<>% as.data.frame()
names(temptable) <- c("inc_med_cat", "freq")
temptable %<>% mutate(
  cumfreq = cumsum(freq)
)
temptable %<>% mutate(
  income_dec = case_when(
    cumfreq <= 0.1 ~ 1,
    cumfreq > 0.1 & cumfreq <= 0.2 ~ 2,
    cumfreq > 0.2 & cumfreq <= 0.3 ~ 3,
    cumfreq > 0.3 & cumfreq <= 0.4 ~ 4,
    cumfreq > 0.4 & cumfreq <= 0.5 ~ 5,
    cumfreq > 0.5 & cumfreq <= 0.6 ~ 6,
    cumfreq > 0.6 & cumfreq <= 0.7 ~ 7,
    cumfreq > 0.7 & cumfreq <= 0.8 ~ 8,
    cumfreq > 0.8 & cumfreq <= 0.9 ~ 9,
    cumfreq > 0.9 ~ 10
  )
)

temp %<>% 
  left_join(temptable[, c("inc_med_cat", "income_dec")], by="inc_med_cat")
temp %<>% group_by(zipcode) %>% slice(1) %>% ungroup()
mediamark %<>% left_join(temp[,c("zipcode", "income_dec")], by="zipcode")





##### Step 12: regression tables #####



mediamark %<>% mutate(
  stateyear = as.factor(paste0(str_sub(as.character(zipcode), 1, -4), "~", as.character(int_year)))
)

# express hh_income in thousands
mediamark %<>% mutate(
  hh_income = hh_income/1000
)

# remove imputed zero hours
mediamark %<>% mutate(
  fnc_hours = ifelse(fnc_hours==0, NA_integer_, fnc_hours)
)

# turn ideology into factor
mediamark %<>% mutate(
  ideology_quintile_1 = (ideology_quintile==1),
  ideology_quintile_0 = (ideology_quintile==0)
)

cable_vars <- 
  c("cnn_hours", "cnn_any", "fnc_hours", "fnc_any", "msnbc_hours", "msnbc_any")
indiv_controls <- 
  c("hh_income", 
    "ideology_quintile",
    "age_quintile",
    "white", "black", "hispanic", "college_degree", "man")
cable_controls <- c(
  "i(chan_config)",
  "num_channels",
  "num_bc")
zip_controls2 <- intersect(zip_controls, names(mediamark))
zip_controls2 <- c(zip_controls2, "i(income_dec)")

for(var in indiv_controls) {
  
  mediamark[[paste0("fnc_pos_interact_", var)]] <- mediamark[["pos_lin_fnc"]]*mediamark[[var]]
  mediamark[[paste0("fnc_hours_interact_", var)]] <- mediamark[["fnc_hours"]]*mediamark[[var]]
  mediamark[[paste0("fnc_any_interact_", var)]] <- mediamark[["fnc_any"]]*mediamark[[var]]
  
  mediamark[[paste0("msnbc_pos_interact_", var)]] <- mediamark[["pos_lin_msnbc"]]*mediamark[[var]]
  mediamark[[paste0("msnbc_hours_interact_", var)]] <- mediamark[["msnbc_hours"]]*mediamark[[var]]
  mediamark[[paste0("msnbc_any_interact_", var)]] <- mediamark[["msnbc_any"]]*mediamark[[var]]
}

indiv_controls_pos_interact <- 
  mediamark %>% dplyr::select(starts_with("fnc_pos_interact")) %>% names()
indiv_controls_hours_interact <- 
  mediamark %>% dplyr::select(starts_with("fnc_hours_interact")) %>% names()
indiv_controls_any_interact <-
  mediamark %>% dplyr::select(starts_with("fnc_any_interact")) %>% names()

indiv_controls_pos_interact_msnbc <- 
  mediamark %>% dplyr::select(starts_with("msnbc_pos_interact")) %>% names()
indiv_controls_hours_interact_msnbc <- 
  mediamark %>% dplyr::select(starts_with("msnbc_hours_interact")) %>% names()
indiv_controls_any_interact_msnbc <-
  mediamark %>% dplyr::select(starts_with("msnbc_any_interact")) %>% names()




for(dv in c("fnc_any", "fnc_hours")) {
  
  print(paste0("-----", dv, "-----"))
  
  f_base <- fixest_form(y=dv, x=c("pos_lin_fnc", "pos_lin_msnbc", cable_controls), fe="0")
  f_controls <- fixest_form(y=dv, x=c("pos_lin_fnc", "pos_lin_msnbc", cable_controls, indiv_controls), fe="0")
  f_controls_zip <- fixest_form(y=dv, x=c("pos_lin_fnc", "pos_lin_msnbc", cable_controls, indiv_controls, zip_controls2), fe="0")
  f_controls_stateyearfe <- fixest_form(y=dv, x=c("pos_lin_fnc", "pos_lin_msnbc", cable_controls, indiv_controls, zip_controls2), fe="stateyear")
  
  
  m_base <- feols(f_base, data=mediamark, cluster="stid")
  m_controls <- feols(f_controls, data=mediamark, cluster="stid")
  m_controls_zip <- feols(f_controls_zip, data=mediamark, cluster="stid")
  m_controls_stateyearfe <- feols(f_controls_stateyearfe, data=mediamark, cluster="stid")
  
  if(dv == "fnc_hours") {
    plottitle = "FNC Channel Positions on Cable Viewership (Intensive Margin)"
  }
  if(dv == "fnc_any") {
    plottitle = "FNC Channel Positions on Cable Viewership (Extensive Margin)"
  }
  
  if(length(m_base$obs_selection$obsRemoved)==0) {
    mean_y1 <- mean(mediamark[[dv]], na.rm=TRUE)
  } else {
    mean_y1 <- mean(mediamark[[dv]][setdiff(1:nrow(mediamark), -m_base$obs_selection$obsRemoved)], na.rm = TRUE)
  }
  mean_y2 <- mean(mediamark[[dv]][setdiff(1:nrow(mediamark), -m_controls$obs_selection$obsRemoved)], na.rm = TRUE)
  mean_y3 <- mean(mediamark[[dv]][setdiff(1:nrow(mediamark), -m_controls_zip$obs_selection$obsRemoved)], na.rm = TRUE)
  mean_y4 <- mean(mediamark[[dv]][setdiff(1:nrow(mediamark), -m_controls_stateyearfe$obs_selection$obsRemoved)], na.rm = TRUE)
  
  setFixest_dict(c(pos_lin_fnc="FNC Channel Pos.",
                   pos_lin_msnbc = "MSNBC Channel Pos.",
                   fnc_hours = "FNC Hours",
                   fnc_any = "Any FNC",
                   stateyear = "State-Year",
                   hh_income = "HH Income (in Thousands)",
                   ideology_quintile = "Conservativism",
                   age_quintile = "Age (Quintile)",
                   white = "White",
                   black = "Black",
                   hispanic = "Hispanic",
                   college_degree = "College Degree",
                   man = "Man"
  ),
  reset=TRUE
  )
  
  setFixest_etable(digits = 3,
                   fitstat = ~ n + r2)
  
  # Print Tables C.1.1 and C.1.2 -----
  print(etable(m_base, m_controls, m_controls_zip, 
         m_controls_stateyearfe,
         extralines = list("__Mean of DV" = c(sprintf("%.3f", mean_y1), sprintf("%.3f", mean_y2),
                                              sprintf("%.3f", mean_y3), sprintf("%.3f", mean_y4))),
         # file = paste0("mediamark-compliance-", dv, ".tex"),
         title = plottitle,
         cluster = ~stid,
         digits=3,
         digits.stats=2,
         drop = "(Constant)",
         label = paste0("table:mediamark-compliance-", dv),
         style.tex=tablestyle,
         group=list("\\midrule Cable System Characteristics"=c("chan", cable_controls),
                    "Zip Code Demographics"=c("income_dec", zip_controls2)),
         replace = TRUE,
         placement="H"
  )) 
  
}



for(dv in c("fnc_any", "fnc_hours")) {
  
  print(paste0("-----", dv, "-----"))
  
  f_base <- fixest_form(y=dv, x=c("pos_lin_fnc", "pos_lin_msnbc", cable_controls), fe="0")
  f_controls <- fixest_form(y=dv, x=c("pos_lin_fnc", "pos_lin_msnbc", indiv_controls,
                                      indiv_controls_pos_interact, indiv_controls_pos_interact_msnbc, 
                                      cable_controls), fe="0")
  f_controls_zip <- fixest_form(y=dv, x=c("pos_lin_fnc", "pos_lin_msnbc", indiv_controls,
                                          indiv_controls_pos_interact, indiv_controls_pos_interact_msnbc, 
                                      cable_controls, zip_controls2), fe="0")
  f_controls_stateyearfe <- fixest_form(y=dv, x=c("pos_lin_fnc", "pos_lin_msnbc", indiv_controls,
                                                  indiv_controls_pos_interact, indiv_controls_pos_interact_msnbc, 
                                             cable_controls, zip_controls2), fe="stateyear")
  
  m_base <- feols(f_base, data=mediamark, cluster="stid")
  m_controls <- feols(f_controls, data=mediamark, cluster="stid")
  m_controls_zip <- feols(f_controls_zip, data=mediamark, cluster="stid")
  m_controls_stateyearfe <- feols(f_controls_stateyearfe, data=mediamark, cluster="stid")
  
  if(dv == "fnc_hours") {
    plottitle = "FNC Channel Positions on Cable Viewership (Intensive Margin)"
  }
  if(dv == "fnc_any") {
    plottitle = "FNC Channel Positions on Cable Viewership (Extensive Margin)"
  }
  
  setFixest_dict(c(
                   pos_lin_fnc="FNC Channel Pos.",
                   pos_lin_msnbc = "MSNBC Channel Pos.",
                   fnc_hours = "FNC Hours",
                   fnc_any = "Any FNC",
                   stateyear = "State-Year",
                   fnc_pos_interact_hh_income = "FNC Channel Pos. $\\times$ HH Income",
                   fnc_pos_interact_ideology_quintile = "FNC Channel Pos. $\\times$ Conservatism",
                   fnc_pos_interact_age_quintile = "FNC Channel Pos. $\\times$ Age",
                   fnc_pos_interact_white = "FNC Channel Pos. $\\times$ White",
                   fnc_pos_interact_black = "FNC Channel Pos. $\\times$ Black",
                   fnc_pos_interact_hispanic = "FNC Channel Pos. $\\times$ Hispanic",
                   fnc_pos_interact_college_degree = "FNC Channel Pos. $\\times$ College Degree",
                   fnc_pos_interact_man = "FNC Channel Pos. $\\times$ Man"
  ),
  reset=TRUE
  )
  
  if(length(m_base$obs_selection$obsRemoved)==0) {
    mean_y1 <- mean(mediamark[[dv]], na.rm=TRUE)
  } else {
    mean_y1 <- mean(mediamark[[dv]][setdiff(1:nrow(mediamark), -m_base$obs_selection$obsRemoved)], na.rm = TRUE)
  }
  mean_y2 <- mean(mediamark[[dv]][setdiff(1:nrow(mediamark), -m_controls$obs_selection$obsRemoved)], na.rm = TRUE)
  mean_y3 <- mean(mediamark[[dv]][setdiff(1:nrow(mediamark), -m_controls_zip$obs_selection$obsRemoved)], na.rm = TRUE)
  mean_y4 <- mean(mediamark[[dv]][setdiff(1:nrow(mediamark), -m_controls_stateyearfe$obs_selection$obsRemoved)], na.rm = TRUE)
  
  # Print Tables C.2.1 and C.2.2 -----
  print(etable(m_base, m_controls, m_controls_zip, 
         m_controls_stateyearfe,
         extralines = list("__Mean of DV" = c(sprintf("%.3f", mean_y1), sprintf("%.3f", mean_y2),
                                              sprintf("%.3f", mean_y3), sprintf("%.3f", mean_y4))),
         # file = paste0("mediamark-compliance-heterogeneity-", dv, ".tex"),
         title = plottitle,
         cluster = ~stid,
         digits=3,
         digits.stats=2,
         drop = "(Constant)",
         label = paste0("table:mediamark-compliance-heterogeneity-", dv),
         style.tex=tablestyle,
         group=list("\\midrule Cable System Characteristics"=c("chan", cable_controls),
                    "Respondent Demographics"=c(indiv_controls, indiv_controls_pos_interact_msnbc),
                    "Zip Code Demographics"= c("income_dec", zip_controls2)),
         replace = T,
         placement="H"
  ) )
  
  
}




